home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Libris Britannia 4
/
science library(b).zip
/
science library(b)
/
PROGRAMM
/
CC_C
/
0924.ZIP
/
SMFPS.LIB
< prev
next >
Wrap
Text File
|
1988-01-10
|
9KB
|
236 lines
/* SMC float function library */
#list-
/* LOGARTIHMIC AND TRIGONOMETRICAL FUNCTIONS */
/* Trig functions are all in degrees measure */
float exp (x) float x;
{ float fint(), fabs(), chebgen();
static float e, m, n, expcc[] =
{ 1.4569998748 E+0, 2.4876243388 E-1, 2.1446555995 E-2,
1.2357140820 E-3, 5.3453058176 E-5, 1.8506907121 E-6,
5.3411876864 E-8, 1.3215163808 E-9, 2.8613237815 E-11 };
static int v, *ex = &e, *mx = &m,
ns = sizeof (expcc) / sizeof (float);
if (fabs (x) >= 400.) return (x < 0.) ? 0. : 4. E+153;
v = n = fint (m = x * 1.4426950406);
m -= n; if (*mx) ++(*mx); e = chebgen (m - 1., expcc, ns);
v = (*ex += v); if (v >= 0x400) return 4. E+153;
return (v <= 0) ? 0. : e; }
float loge (x) float x;
{ float chebgen();
static float e, m, n, logcc[] =
{ 9.3022922134 E-1, -8.1841456644 E-2, 9.4766115856 E-3,
-1.2282837393 E-3, 1.6939529027 E-4, -2.4301272033 E-5,
3.5827616763 E-6, -5.3890849189 E-7, 8.2315218987 E-8,
-1.2726716189 E-8, 1.9871237790 E-9, -3.1280003518 E-10,
4.9576937640 E-11 };
static int v, *nx = &n, ns = sizeof (logcc) / sizeof (float);
if (x <= 0.) return -4. E+153;
n = x; v = *nx - 0x200; *nx = 0x200;
if (n <= 0.8) { --v; ++(*nx); }
e = chebgen (n * 2.5 - 3.0, logcc, ns);
return e * (n - 1.) + (v * 0.69314718025); }
float tan (x) float x;
{ float m, _stgen(), sqrt();
m = _stgen (x, 1); return m / sqrt (1. - m * m); }
float cos (x) float x;
{ float _stgen(); return _stgen (x + 90., 0); }
float sin (x) float x;
{ float _stgen(); return _stgen (x, 0); }
float _stgen (x, t) float x;
{ float fint(), chebgen();
static float m, n, ccsin [] =
{ 1.2762789622 E+0, -1.4263078457 E-1, 4.5590080009 E-3,
-6.8293756748 E-5, 5.9248092843 E-7, -3.3513958000 E-9,
1.3336392981 E-11 };
static int *nx = &n, ns = sizeof (ccsin) / sizeof (float);
n = x/360.; n -= fint (n + 0.5); if (*nx) *nx += 2;
if ((*nx & 0x7FFF) > 0x200) n = t ?
(n > 0. ? -2. : 2.) + n : (n > 0. ? 2. : -2.) - n;
n *= (m = n); if (*nx) ++(*nx); --n;
return m * chebgen (n, ccsin, ns); }
float acos (x) float x;
{ float atan(), sqrt();
return 90. - atan (x / sqrt (1. - x * x)); }
float asin (x) float x;
{ float atan(), sqrt();
return atan (x / sqrt (1. - x * x)); }
float atan (x) float x;
{ float fabs(), chebgen();
static float m, n, atancc[] =
{ 8.8137358699 E-1, -5.2946462263 E-2, 5.5679210277 E-3,
-6.9059750155 E-4, 9.2871486598 E-5, -1.3107598052 E-5,
1.9105182955 E-6, -2.8495930820 E-7, 4.3244389300 E-8,
-6.6516919840 E-9, 1.0342528811 E-9, -1.6224319589 E-10,
2.5639831258 E-11 };
static int t, *nx = &n, ns = sizeof (atancc) / sizeof (float);
m = (t = (fabs (x) >= 1.)) ? -1./x : x;
n = m * m; if (*nx) ++(*nx);
n = m * 57.295779513 * chebgen (n - 1., atancc, ns);
if (t) n = (x > 0.) ? n + 90. : n - 90.; return n; }
float chebgen (zz, cc, n) float zz, cc[];
{ static float a0, a1, a2, *cp; static int *ax = &a0;
cp = cc + n; a0 = a1 = 0.;
while (cp > cc) { a2 = a1; a1 = a0; a0 = zz * a1;
if (*ax) ++(*ax); a0 += (*--cp - a2); }
return a0 - a2; }
/* GENERAL FLOAT FUNCTIONS */
float fabs (b) float b;
/* return absolute value */
{ static unsigned *e; e = &b;
*e &= 0x7FFF; return b; }
int fsgn (b) float b;
/* return signum: -1, 0 or 1 */
{ static unsigned *e; e = &b;
if (*e) return (*e & 0x8000) ? -1 : 1; return 0; }
float sqrt (x) float x;
/* return square root */
{ if (x <= 0.) return 0.; &x; inline (
0x89, 0xDF, 0x8B, 0x55, 0x04, 0x8B, 0x4D, 0x02, 0x8B, 0x1D, 0x31,
0xF6, 0x89, 0xD8, 0x81, 0xEB, 0xFF, 0x01, 0xD1, 0xFB, 0x81, 0xC3,
0x00, 0x02, 0x25, 0x01, 0x00, 0x74, 0x06, 0xD1, 0xE9, 0xD1, 0xDA,
0xD1, 0xDE, 0x56, 0x52, 0x51, 0x31, 0xC9, 0x31, 0xD2, 0xB8, 0x00,
0x40, 0x51, 0x51, 0x50, 0x89, 0xE5, 0xB8, 0x21, 0x00, 0x50, 0xD1,
0xD0, 0x35, 0x01, 0x00, 0x0B, 0x76, 0x04, 0x0B, 0x56, 0x02, 0x0B,
0x4E, 0x00, 0x29, 0x76, 0x0A, 0x19, 0x56, 0x08, 0x19, 0x4E, 0x06,
0x73, 0x02, 0xD1, 0xD8, 0x73, 0x14, 0x01, 0x76, 0x0A, 0x11, 0x56,
0x08, 0x11, 0x4E, 0x06, 0x33, 0x76, 0x04, 0x33, 0x56, 0x02, 0x33,
0x4E, 0x00, 0xEB, 0x09, 0x03, 0x76, 0x04, 0x13, 0x56, 0x02, 0x13,
0x4E, 0x00, 0xD1, 0x6E, 0x00, 0xD1, 0x5E, 0x02, 0xD1, 0x5E, 0x04,
0xD1, 0x66, 0x0A, 0xD1, 0x56, 0x08, 0xD1, 0x56, 0x06, 0x58, 0x48,
0x75, 0xAF, 0x83, 0xC4, 0x0C, 0x01, 0xF6, 0x73, 0x04, 0x42, 0x75,
0x01, 0x41, 0x89, 0x55, 0x04, 0x89, 0x4D, 0x02, 0x89, 0x1D ); }
float fint (b) float b;
/* return next int value below - eg -1.5 gives -2. */
{ static float a; static unsigned aex;
a = b; aex = * (unsigned *) &a;
if ((aex & 0x7FFF) <= 0x200)
return (aex & 0x8000) ? -1.0 : 0.0;
if ((aex & 0x7FFF) >= 0x220 ) return a;
/* shift down mantissa by (0x220 - exp)
inc if negative and renormalise */
inline ( 0xBD, &a, 0x31, 0xDB, 0xB9, 0x20, 0x02,
0x8B, 0x46, 0x00, 0x25, 0xFF, 0x7F, 0x29, 0xC1,
0xD1, 0x6E, 0x02, 0xD1, 0x5E, 0x04,
0x83, 0xD3, 0x00, 0xFF, 0x46, 0x00,
0x49, 0x75, 0xF1, 0x09, 0xDB, 0x74, 0x11,
0x8B, 0x5E, 0x00, 0x81, 0xE3, 0x00, 0x80, 0x74, 0x08,
0xFF, 0x46, 0x04, 0x75, 0x03, 0xFF, 0x46, 0x02,
0x8B, 0x46, 0x02, 0x0B, 0x46, 0x04,
0x75, 0x05, 0x89, 0x46, 0x00, 0xEB, 0x12,
0xF7, 0x46, 0x02, 0x00, 0x80, 0x75, 0x0B, 0xD1, 0x66,
0x04, 0xD1, 0x56, 0x02, 0xFF, 0x4E, 0x00, 0xEB, 0xEE );
return a; }
/* Convert string representation to float */
float atof (s) char *s;
{ static int e, m; float f, g, fe10();
m = e = 0; f = 0.; g = 1.;
while (*s == ' ' || *s == '\t') ++s;
if (*s == '-' || *s == '+') m = (*s++ == '-');
if (*s && *s != '.')
while (*s >= '0' && *s <= '9')
f = f * 10. + (*s++ - '0');
if (*s == '.') { ++s;
while (*s >= '0' && *s <= '9')
f += (g *= .1) * (*s++ - '0'); }
if (m) f = -f; m = 0;
while (*s == ' ' || *s == '\t') ++s;
if (*s == 'E' || *s == 'e')
{ ++s; while (*s == ' ' || *s == '\t') ++s;
if (*s == '-' || *s == '+') m = (*s++ == '-');
while (*s >= '0' && *s <= '9')
e = e * 10 + (*s++ - '0');
if (m) e = -e; f = fe10 (f, e); }
return f; }
/* PRINTF/SPRINTF/FPRINTF ROUTINES - for FP print capture */
printf (fs)int *fs; {int _psc(), _fprt(); return _prout(stdout, fs, _psc); }
sprintf(fs)int *fs; {int _pss(), _fprt(); return _prout(fs, --fs, _pss); }
fprintf(fs)int *fs; {int putc(), _fprt(); return _prout(*fs, --fs, putc); }
/* Floating point routine for _prout() */
_fprt (f, b, p, t) float f; char *b;
{ if (t) p *= t; else p = (f < 1000000.) ? p : -p;
if ((p = ftoa (f, b, p)) < 0 || t) return;
b += p; while (*--b == '0') *b = ' '; }
/* FP to Ascii representation. Exponential format used
if precision is negative, or abs (f) exceeds 1. E 9.
Returns +/- string length, guaranteed less than 30 chars. */
ftoa (f, s, p) float f; char *s;
{ static int dd[3], ee, ef, sg, zz, qq;
static char *t; float fe10();
struct { unsigned _exp, _mhi, _mlo; };
if ((qq = p) < 0) { ef = 1; p = -p; } else ef = 0;
if (p > 14) p = 14; sg = f._exp & 0x8000;
if ((f._exp &= 0x7FFF) > 0x21E) ef = 1;
if (f._exp)
{ if (!ef) f += fe10 (.5, -p);
ee = f._exp - 0x202;
ee = (3 * ee + ee / 100) / 10 + (ee > 0 ? 1 : 0);
f = fe10 (f, -ee);
if (f._exp >= 0x201) { f = f * .1; ++ee; }
if (ef)
{ if (qq >= 0) p += 2; f += fe10 (.05, -p);
if (f._exp >= 0x201) { f = f * .1; ++ee; } }
dd[0] = f._mhi; dd[1] = f._mlo; dd[2] = 0;
if ((zz = 0x200 - f._exp) > 0)
while (zz--) /* shift array right */
inline ( 0xBB, dd, 0xD1, 0x2F, 0x43, 0x43,
0xD1, 0x1F, 0x43, 0x43, 0xD1, 0x1F ); }
else { ee = 1; dd[0] = dd[1] = dd[2] = 0; }
t = s + 1; zz = 0;
if (!ef) { if (ee <= 0) *t++ = '0';
else do *t++ = '0' + _fpc (dd, &zz); while (--ee); }
else *t++ = '0' + _fpc (dd, &zz); *t++ = '.';
if (!ef && ee) while (ee++ && p-- > 0) *t++ = '0';
while (p-- > 0) *t++ = '0' + _fpc (dd, &zz);
if (ef) { *t++ = ' '; *t++ = 'E';
if (--ee >= 0) *t++ = '+';
else { *t++ = '-'; ee = -ee; }
if (ee >= 100) { *t++ = '0' + ee/100; ee %= 100; }
*t++ = '0' + ee/10; *t++ = '0' + ee % 10; }
*s = (zz && sg) ? '-' : ' '; *t = 0;
return ef ? s - t : t - s; }
/* Unload characters from array - or with c */
_fpc (a, c) int *a, *c;
{ inline ( 0x89, 0xE5, 0x8B, 0x76, 0x04, 0x83, 0xC6, 0x06,
0x31, 0xDB, 0xB9, 0x03, 0x00, 0x4E, 0x4E, 0xB8,
0x0A, 0x00, 0xF7, 0x24, 0x01, 0xD8, 0x89, 0x04,
0xBB, 0x00, 0x00, 0x11, 0xD3, 0xE2, 0xEE,
0x8B, 0x76, 0x02, 0x09, 0x1C ); }
float fe10 (b, p) float b;
/* return b * 1. E+p */
{ static float m;
if (p > 0) m = 10.; else { p = -p; m = .1; }
while (p) { if (p & 0x1) b *= m;
if (p >>= 1) m *= m; } return b; }
/* End of float library */
#list+